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ABSTRACT 

Mercury's eccentricity is chaotic and can increase so much that collisions with Venus or the Sun become possible. This chaotic 
behavior results from an intricate network of secular resonances, but in this paper, we show that a simple integrable model with 
only one degree of freedom is actually able to reproduce the large variations in Mercury's eccentricity, with the correct amplitude 
and timescale. We show that this behavior occurs in the vicinity of the separatrices of the resonance gi - g 5 between the precession 
frequencies of Mercury and Jupiter. However, the main contribution does not come from the direct interaction between these two 
planets. It is due to the excitation of Venus' orbit at Jupiter's precession frequency g s . We use a multipolar model that is not expanded 
with respect to Mercury's eccentricity, but because of the proximity of Mercury and Venus, the Hamiltonian is expanded up to order 
20 and more in the ratio of semimajor axis. When the effects of Venus' inclination are added, the system becomes nonintegrable and 
a chaotic zone appears in the vicinity of the separatrices. In that case, Mercury's eccentricity can chaotically switch between two 
regimes characterized by either low-amplitude circulations or high-amplitude librations. 

Key words, methods: analytical - methods: numerical - chaos - celestial mechanics - planetary systems - planets and satellites: 
dynamical evolution and stability - planets and satellites: individual: Mercury 



1. Introduction 

After the discovery of the chaotic motion of the planets in the 
Solar System (Laskar 1989, 1990), it has been demonstrated 
that the eccentricity of Mercury can rise to very high values 
(Laskar 1994), allowing for collision of the planet with Venus. 
This was confirmed later on by direct numerical integration in 
simplified models that neglect general relativity (GR) (Laskar 
2008; Batygin & Laughlin 2008), and even in the full model 
that includes the GR contribution (Laskar & Gastineau 2009). 
Quite surprisingly, the behavior of the system depends strongly 
on the GR contribution. Indeed, the probability that the eccen- 
tricity increases beyond 0.7 in less than 5 Gyr is about 1% in 
the full model, while it rises to more than 60% when GR is ne- 
glected (Laskar 2008; Laskar & Gastineau 2009). This behavior 
is due to the presence of a secular resonance between the peri- 
helion motions of Mercury (vji) and Jupiter {vjs) (Laskar 2008; 
Batygin & Laughlin 2008). Although the origin of the instability 
in Mercury's eccentricity is identified, the precise mechanism of 
the behavior of Mercury's orbit has not yet been explained in 
detail. A first attempt was provided by (Lithwick & Wu 2011), 
who analyzed the overlap of resonances in the truncated secu- 
lar Hamiltonian of degree four in eccentricity and inclination. 
This model reproduces the small diffusion of Mercury's eccen- 
tricity well, and confirms the important role of the gi — gs and 
(gi _ gs) _ ( s i - S2) secular resonances between the precession 
frequency of the perihelion (g,) and the regression rate of the 
ascending node (5,). However, Lithwick & Wu (2011) does not 
explain the steady increase in Mercury's eccentricity beyond 0.7 
as observed in (Laskar 2008). 

In the present paper, our approach is different, since we do 
not want to be bound to the limitations given by expansions 



in eccentricity as in (Lithwick & Wu 2011) or in the previous 
work of Laskar (1984, 1989, 1990, 2008). As we are looking 
for a large excursion of the eccentricity of Mercury, we prefer 
to derive expressions that are valid for all eccentricities, using 
the averaged expansions derived in (Laskar & Boue 2010). We 
provide two simple analytical, secular models. The first one is 
coplanar with only one degree of freedom. It is thus integrable. 
This model shows that Mercury's eccentricity can reach values 
as high as 0.8 if the system is in the vicinity of the g\ - gs res- 
onance. The second model includes inclination and has two de- 
grees of freedom. This spatial model exhibits the chaotic behav- 
ior of Mercury's eccentricity in the neighborhood of the separa- 
trix with possibilities of switching between a regime of moderate 
eccentricity to a regime of large oscillation where Mercury's ec- 
centricity reaches very high values. In both models, Mercury 
is treated as a massless particle, while the motion of the other 
planets in the Solar System are given by a single term in their 
quasiperiodic decompositions taken from (Laskar 1990). In sim- 
plified terms, the chaos in Mercury's orbit is due to the proximity 
of the resonance with the precession motion of Venus excited at 
the frequency of Jupiter's precession frequency. 

The secular models that we use in this paper rely on a multi- 
polar expansion of the perturbing function, up to order (ai/a p ) 20 
in the relativistic case and (ai/a p ) 50 in the Newtonian case, 
where (a p ) p= 2,& are the semimajor axes of the planets of the Solar 
System in increasing order. This expansion, subsequently aver- 
aged over the mean anomalies of all the planets, allows for arbi- 
trary inclinations and eccentricities as long as no orbit crossing 
occurs. 

In section 2, we derive the equations of motion of the copla- 
nar model. This model is very simple with only one degree of 
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freedom associated to Mercury's eccentricity and longitude of 
perihelion. Then, in section 3, we derive the possible trajecto- 
ries in the phase space using level curves of the Hamiltonian 
for different orders of expansion. We show that it is necessary 
to develop the perturbing function up to high orders in (a\la p ) 
in order to reach the asymptotic evolution. We also show that 
the maximum eccentricity attained with an initial eccentricity 
of ei = 0.2 is on the order of 0.8 within 4 Myrs as reported 
in (Laskar 2008). The spatial case is treated in section 4. This 
nonintegrable model illustrates the generation of chaos in the 
vicinitiy of the hyperbolic fixed point of the planar model. We 
conclude in the last section. 



2. Coplanar model 

2.1. Newtonian interaction 

Using the traditional notations where planets have increasing in- 
dices with respect to their semimajor axis, we note the barycen- 
tric position of the Sun Ho, of Mercury U\, of Venus 112, etc. The 
heliocentric positions of the planets are similarly noted (r p ) p= i 
Since the mass of Mercury is much lower than the masses of 
the other planets in the Solar System, we model Mercury as a 
massless particle. 

Our goal is to study the behavior of Mercury under a known 
planetary perturbation. As such, all of the u p and their deriva- 
tives, for p + 1, are considered as given functions of the time 
t. Using the Poincare heliocentric canonical variables (Laskar 
& Robutel 1995), the Hamiltonian describing the evolution of 
Mercury's trajectory reads as 1 
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— h-Up(t) 
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(1) 



where f\ = Mi is the conjugate momentum of r\, G the gravita- 
tional constant, mo the mass of the Sun, and {m p ) p= 2$ the masses 
of the other planets. The first two terms of this Hamiltonian rep- 
resent the Keplerian motion and are equal to -Gmo/(2ai), where 
a\ is the semimajor axis of Mercury. 

Considering that Mercury is the closest planet to the Sun, we 
now expand formally as a series in (ri/r p ) p= 2,s 



Hk — — 



Gmo 
2a i 



-cZ?Z 



=2 'P n=0 \'P 
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where P„ is the Legendre polynomial of order n. 

A first approximating model could stop the expansion at the 
second Legendre polynomial Pi- However, it is well known that 
the resulting double-averaged secular Hamiltonian does not de- 
pend on the perihelia of the outer bodies at this order (Lidov 
1962;Kozai 1962; Lidov &Ziglin 1976; Farago & Laskar 2010), 
so no secular resonance between the perihelia of Mercury and 
any other planet can be seen there. We thus push the expansion 
to higher orders. Furthermore, we see in the following that it is 
necessary to extend the sum at least up to n ~ 10. 

The next step consists in averaging the Hamiltonian (2) over 
the mean anomalies (Mp)p=i,s °f all planets. In the Hamiltonian 
(2), the terms —Gmo/r p obtained for n — do not depend on 
Mercury's elements; the terms (f\ ■ u p ) vanish when averaged 
over M\, the term — Gmo/(2a 1 ) becomes constant after averag- 
ing over Mi, since a\ becomes constant. As such, all these terms 



1 For clarity, we drop the explicit time dependence in all equations 
after (1). 



are dropped in the following expressions. Finally, the averaged 
expression of the perturbing functions a p / \r\ - rJ are given in 
(Laskar & Boue 2010). The resulting Hamiltonian of the copla- 
nar problem is then 
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In these expressions, a p , e p , and m p are the semimajor axis, the 
eccentricity, and the longitude of the pericenter of the planet p, 
respectively. The X'^ m (e) are the Hansen coefficients defined by 



oo 

-j e =}_ j X k (e) e , 



&=-oo 



where e„ — if n is odd, and e„ — 1 if n is even, and 
(2q)\(2n - 2q)\ 



fn,q — 
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To simplify the expansion of the Hamiltonian (2), we take 
advantage of the eccentricities {e p ) p= 2,% of all the planets beyond 
Mercury remaining low to only keep the linear terms in these 
eccentricities. In contrast, the expansion remains exact in the ec- 
centricity e\. With this approximation, up to the octupole order, 
the averaged Hamiltonian (2) reads as 



G V~1 ( ®1 2 

p=2 \ a P 



(7) 
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More generally, as long as the Hamiltonian is truncated at 
the first order in e p , its expression can be written as (see 
Appendix A) 



H, 
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where PJa, x) and Q,^a, x) are polynomials of degree [n/2] and 
[(n - l)/2] in x and of degree 2[n/2] and 2[(n + l)/2] - 1 
in a, with n the degree of expansion of the perturbing func- 
tions as in Eq. (3). It can be noted that the coefficient of x° in 
Po(a, x) is the Taylor expansion of C\{a) - 1 = (l/2)by 2 (a) 
-1, and the coefficient of x 1 in PJ^a, x) is the Taylor series of 
Ci(a)l2 = (l/S)ab^ 2 (a), etc., where bf\a) are Laplace coeffi- 
cients (Laskar & Robutel 1995). 

In the Hamiltonian (8), e p and m p are functions of time. 
Neglecting the slow diffusion of the eccentricities of the inner 
planets, the evolution of each z p = e p exp im p is described well 
by a quasiperiodic expansion (Laskar 1990). In this study, we 
focus on the secular resonance g t — g$. In the quasiperiodic de- 
composition of the variables (z p ) p =2,$, we thus keep only the 
terms associated to the eigenmode z* = exp igst with frequency 
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ei cos A-ca e\ cos Aro cos Aro 



Fig. 1. Level curves of the Hamiltonian in the plane {e\ cos Aw, e\ sin Act), a) Quadrupole expansion with g\ = 5"/yr. b) Octupole 
expansion with gi = 5"/yr. c,d,e,f) expansion up to the order n — 20 with g\ = 2.5, 3.6, 4.0, 5.0"/yr, respectively. 



Table 1. Amplitude of the quadrupole £^ uad (14) and octupole 
g octu 3) terms, due to each planet p, amplitude A p , and phase 
<p p of the eigenmode with frequency gs = 4.2488 163"/yr in the 
quasiperiodic decomposition of z p - e p exp im p (taken from 
Laskarl990). 
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4 uad x io 6 


e octu x 1Q 6 


A„ x 10 6 


<P P (deg.) 
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59 375 


1 170 


19 636 


30.571 


3 


27911 


383 


18913 


30.597 


4 


837 
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20300 


30.679 


5 


62201 


383 


44119 


30.676 


6 


3 025 


8 


33 142 


30.676 


7 


57 





37 351 


210.671 


8 


17 





1771 


30.669 



Gm p /(Aia p ) = n\a\m p j(a p m^) , the resonant Hamiltonian now 
reads as 



m„ ail I a \ i\ 

H N , p i dn = -«i > Pq. —,e l 

j-^ mi) a p \ \a p ) 



-ei A p Q a \ — , e? I cos(cr 1 - vj*) 



(9) 



With this convention, the time should also be rescaled by the 
same factor Ai to keep the canonical form of the equations of 
motion, unless the canonical variables are modified as follows. 
We let 



7 = A 



(10) 



be the canonical conjugated variables of the initial Hamiltonian 
#jv,pian °f (Eq.8). It can be easily shown that 



gs ~ 4.25"/yr. The others average out and disappear from the 
Hamiltonian. 

The amplitude and the phase of the mode z* in the decom- 
position of each variable z p are provided in Table 1 . We notice 
that, except for Uranus (p = 7), all the phases are the same and 
equal to « 30.6 deg. Since, tp-i = 30.6 + 180 deg, it is equiva- 
lent to take ifL = 30.6 deg and A' 7 = —Aj, This is the convention 
that is followed hereafter, but the prime is omitted for clarity. 
We note zu* = gst + <p*, where ip* = 30.6 deg. We also use the 
fact that A! = -\JGmoO] is constant in the secular problem to 
rescale the Hamiltonian by this quantity as it simplifies the fol- 
lowing expressions. We note Z/jv.pian = ffjv.pianM-i • Given that 
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are canonical conjugated variables of the Hamiltonian ff;v,pian 
(Eq.9) without rescaling the time t. The equations of motion 
are thus 



dl 
dt 



d J^iiXt), 

88 



d6 dH Np i m y v. 

-r = — e—(i,0,t) 

dt 81 



(12) 



At this point, it is interesting to evaluate the contribution of 
each planet to the octupole interaction with Mercury. Although 
the resonance g\ - gs involves only the precession frequencies 
associated to Mercury and Jupiter, the eigenmode z* is present 
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Fig. 2. Evolution of the fixed points as a function of the precession frequency of Mercury's orbit at zero eccentricity and for 
different expansion orders. Dashed and dotted curves correspond to unstable points, while the solid ones show the positions of 
the stable points. The equilibrium points are labeled with roman letters as in Fig. 1. The thin vertical line indicates the resonance 

gl = 85- 



in the quasiperiodic decomposition of the motion of all the plan- degree of Q(a, x) in a is 3, the contribution of the octupole terms 
ets. Furthermore, the amplitudes of this mode are very similar can be estimated using the parameters (e°/ tu ) p =2.8 given by 
and vary only within a factor 2.5 between 0.019 and 0.044, ex- 
cept for Neptune whose amplitude is only 0.002 (Table 1). From oc(u 15 ln\\lm p \la\ \ 4 
the expression of the Hamiltonian (9), knowing that the lowest 

£ ° ,cu= 64 y % w Ap ' (13) 
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The values taken by these parameters are gathered in Table 1. 
The maximal amplitude is due to Venus with e^ 1 " * 0.0012 fol- 
lowed by the Earth-Moon barycenter and Jupiter with e° ctu * 
e octu ~ 4 x jo -4 . Thus, the strongest perturbation on Mercury's 
orbit comes from the precession of Venus excited by Jupiter. 
ai/fl2 ~ 0.53 is not very small explains why it is necessary to 
perform the expansion of the perturbing function up to a high or- 
der in the ratio of the semimajor axes. For completeness, Table 1 
also provides the quadrupole contribution of each planet given 
by 



quad 



1 In 



g5/\mo 



(14) 



Since we have shown that several planets play a role in the evo- 
lution of Mercury's eccentricity, in the following we always con- 
sider all the perturbers from Venus to Neptune. To simplify the 
notations, we define two new polynomials P and Q of the eccen- 
tricity e\ alone, given by (Eq. 9) 
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p=2 
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Qw = - y 

g 5 m a p 



m Q a p 
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Then, the resonant Hamiltonian (9) reads as 



77, 



/V,plan 



-f S (P(*l) - «lQ(«?) COSC^J - VT*)) 



(16) 



The numerical values of the coefficients of the polynomials P 
and Q are given in Appendix B. 

The Hamiltonian 7?/v,pian (16) has one and a half degrees 
of freedom, but it can be reduced to a one degree of freedom 
Hamiltonian fljv.pian after some modifications. For that, we first 
make the Hamiltonian autonomous by a adding a momentum T 
conjugated to time t. Then, we perform a canonical transforma- 
tion where the old variables are 



(/=i- ^i-e\,e = -m); (f,t), 

and the new ones are (7, Am) and (f, t), defined by 
9-g 5 t-(p* = m\ -m* , 



Am ■ 
t= t 



(17) 



(18) 



For this transformation to be canonical, the momentums should 
verify 



/=-/=i- a/i-' 

f = -g 5 I + f. 



(19) 



The new Hamiltonian 77^ p i an expressed in the new variables does 
not depend on the cyclic coordinate I. Thus, its conjugated mo- 
mentum T is an integral of the motion and can be dropped. Up 
to a constant, the new Hamiltonian is then 



#/v,pi a n = ~gs yl ~e\ - g 5 P (el) + g5«iQ (efj cos Am . 



(20) 



This one degree of freedom Hamiltonian is integrable. The or- 
bits are given by Hji/^im - Cte. The temporal evolution of the 



eccentricity e\ and of the resonant angle Am are deduced from 
the canonical equations of motion. This leads to 



de\ 
dt 

dAm 
dt 
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2.2. General relativistic precession 

The secular effect of relativity is described by (e.g., Touma et al. 
2009), 



Hr = -g r 



1 



(22) 



where g r = 3(Gmo) 2 /(AiOjC 2 ), and c is the speed of light. In the 
case of Mercury, we have g r as 0.41"/yr. The total Hamiltonian 
^R.pian = ffiv.pian + H R , including the Newtonian interaction and 
general relativity, becomes 



77, 
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-gsl + P (4) ~ eiQ (ef) cos Am) 



(23) 



5r / ' 



2.3. Additional control term 

Using the Newtonian Hamiltonian (20), or the relativistic one 
(23), the system is far from the g\ - gs resonance, in a con- 
figuration where the amplitude of oscillation of the eccentricity 
is small (see next section). Indeed, taking an order of expan- 
sion n = 50, we get g\ = 5.54"/yr in the Newtonian case and 
gi = 5.96"/yr in the relativistic case, whereas the resonance oc- 
curs in the vicinity of gi — g 5 = 4.25"/yr. This result is slightly 
different, but representative of the present behavior of Mercury's 
eccentricity. Indeed, the present value of g\ for the Solar System 
with GR is gi = 5.59"/yr (Laskar et al. 2004), the differences 
with the present model being due to the simplifications that are 
made here. To recover a model that is dynamically close to 
the real Solar System, we add a correction to the Hamiltonian 
which changes the value of gi by an increment 6 g . In addition, 
owing to the slow chaotic diffusion in the inner Solar System, 
Mercury's precession frequency g t can change quasi-randomly 
and come close to the value of gs, then leading to high unsta- 
bility (Laskar 1990, 2008; Batygin & Laughlin 2008; Laskar & 
Gastineau 2009). To explore the evolution of Mercury's eccen- 
tricity for different values of g\ around the resonant frequency 
g$, the above-mentioned correction is added to both ffjv.pian an d 

#fi,plan; i-C, 



-Kaplan 

and 



77, 



1 



N.plan 



#fi,plan - #fi,plan + 2^S e l 



(24) 



(25) 



In Eqs. (24) and (25), the factor 6 g controls the precession fre- 
quency gi. For 6 g = 0, we recover the values obtained with 
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ff/v,pian an d ^R.pian, respectively. Otherwise the value of g\ is in- 
cremented by 5 g . Appendix C provides the values of 5 g that put 
the system in exact resonance for all orders of expansion of the 
Hamiltonians. 



3. Eccentricity behavior 

Here, we analyze the possible trajectories and evolutions of 
Mercury's eccentricity described by the Hamiltonians //jv iP ian 
(24) and ^ (25) obtained in the previous section. Since these 
Hamiltonians are integrable (they have only one degree of free- 
dom), all the orbits are necessarily regular. Thus, the goal of this 
section is not to reproduce the chaotic behavior of Mercury's ec- 
centricity e\ , but to show that e\ can actually reach values as high 
as observed by Laskar (2008) (see Fig. 3a) when the system is in 
the vicinity of the resonance g\ - gs- The chaotic evolution will 
be analyzed in the spatial case (see Sect. 4), where a new degree 
of freedom is added. 



3.1. Dependency with the order of expansion 

The Hamiltonians of the previous section have been computed 
as series in power of the semimajor axis ratios (ai/a p ) p =2,8 for 
any order n. Here, we study the effect of the truncation of these 
series. Hereafter, we note with a superscript (n) any Hamiltonian 
expanded up to the order n, e.g., tr^ l lan or fl^'l . In a first 
step, we focus only on the relativistic Hamiltonian which is more 
complete. 

Within the quadrupole approximation, the Hamiltonian 
#R 2 pian reduces to a simple expression 



H' 



(2) 

i?,plan 



-gs 



,,quad 



(2 + 3efjj 



(26) 



+ 2<Vi 



where e quad = Yjp £ 7> Ud • This Hamiltonian is independent of 
Aw, thus e\ remains constant. All the trajectories in the plane 
(x = e\ cos Am, y — e\ sin Am) are circles centered on the origin 
(0,0) with an eventually infinite period when Am = 0. Figure la 
illustrates this result for gi = 5"/yr. The position of the fixed 
points is represented by red dashed curves. Within this approxi- 
mation, as noted before, there is no possible resonance between 
mi and m* and thus, no possible increase in the eccentricity. 

The next level of approximation is the octupole. The corre- 
sponding Hamiltonian reads as 



H 



(3) 

i?,plan 



\-e\ + e quad (2 + 3e 2 ) 



-e octu (4 + 3ei)cosAzrr|-g,. 



+ <Vi 



(27) 



with e oc,u = Yj P s° ctu . In that case the phase space can be much 
more complicated with five fixed points and two sets of separa- 
trices (see Fig. lb obtained with gi = 5"/yr). The stable fixed 
points are represented by black dots, and the unstable ones are 
noted with crosses. The dashed curves correspond to the separa- 
trices. The positions and existence of the fixed points depend on 
the value of the precession frequency g\. This is depicted in the 
subpanel labeled a 3 of Fig. 2. The curves represent the position 
on the x axis of the fixed points of the phase space as a function 




2 3 4 

time (Myr) 

Fig. 3. Comparison between the eccentricity evolution observed 
in (Laskar 2008) (a), and produced by the simple model of this 
work (b). 



of gi. The labels in roman numerals qualifying each fixed point 
are identical to those in Fig. 1. 

It is interesting to have a look at all of Fig. 2. Indeed, as the 
order n of the expansion of the Hamiltonian H R , increases, 

the degrees of the polynomials P^a, e 2 ) and Qj^a, e 2 ) increase 
both for a and e 2 . Then, the topology of the phase space of the 
Hamiltonian H „ , evolves as shown in Fig. 2. Up to n = 5. 
five fix points - three stable points and two unstable points - 
coexist within 4.5 < g\ < 5.0"/yr, while the system contains at 
most three fix points - two stable points and one unstable point 
- for n > 6. Moreover, the hyperbolic point (labeled IV) located 
at e\ cos Am > 0, i.e. Am = 0, disappears when n > 8. The 
asymptotic topology is reached at n as 10 from which positions 
of the fixed points do not evolve significantly up to n — 20. 

The figures lc, Id, le, and If display the allowed trajectories 
of Mercury's eccentricity for n — 20 which is assumed to be rep- 
resentative of the asymptotic behavior. The values of g\ are 2.5, 
3.6, 4.0, and 5.0"/yr, respectively. At this order n = 20, orbits 
are very similar to those of a simple pendulum with a large res- 
onant island allowing for large increases in eccentricity beyond 
e\ = 0.8. The separatrix disappears at gi > 3.84"/yr but high ec- 
centricity excursions are still possible, to a smaller extent, since 
the elliptic fixed point (V) remains significantly offset with re- 
spect to the center of the phase space as long as g\ < 4.5"/yr. 
For example, with g\ = 4.0"/yr (Fig. le), a trajectory starting 
within <?! < 0.2 can reach e\ at 0.7. 

To conclude, the asymptotic phase space, which must repro- 
duce Mercury's eccentricity behavior as faithfully as possible, 
is obtained for n w 20. Moreover, at this order of expansion, 
Mercury's eccentricity is able to increase from e\ at 0.2 up to 
e\ ai 0.8 as observed in the numerical simulations reported by 
Laskar (2008) (see Fig. 3). 
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Fig. 4. Level curves of the Hamiltonian without relativity in the 
plane (e\ cos Act, e\ sin Am) for two orders of expansion: n = 20 
(left), n = 50 (right). In both cases, g\ = 3.6"/yr. 



3.2. Temporal evolution 



The level curves of the Hamiltonian n£ lw provide the trajecto- 
ries of Mercury's eccentricity in the phase space. In the previous 
section, we saw that the eccentricity is able to vary between 0.2 
and 0.8 as observed in (Laskar 2008) (Fig. 3a). We now check 
whether the timescale of the evolution matches the 4 Myr found 
by Laskar (2008). For this purpose, we integrate the following 
equations 



dx 
dt 

dy 
dt 



- V 



l-e 



dHfl 

2 K,plan 



dy 



■dH 



(20) 
R,plan 



(28) 



dx 



where x — e\ cos Act and y — e\ sin Aw. They are equivalent 
to the equations of motion (21), but without singularities at the 
origin e\ — 0. 

An example of temporal evolution given by the numerical 
integration of Eq. (28) is displayed in Fig. 3b. The initial con- 
ditions are e\ = 0.2, Act = HOdeg, and gi = 3.68"/yr. The 
comparison of Fig. 3a with Fig. 3b shows that the simple one 
degree of freedom model described in this study is able to repro- 
duce the main evolution of Mercury's eccentricity with the cor- 
rect amplitude and timescale. Only once the eccentricity reaches 
a value close to 0.8, Mercury's orbit becomes highly unstable 
due to close encounters with Venus, and the evolution observed 
in (Laskar 2008) starts to differ significantly from the simple 
model. It should be noted that the small oscillations present in 
the integration of the full Solar System (Fig. 3a) do not appear in 
the simple model since we keep only one single eigen frequency 
to describe the evolution of the outer planets eccentricities. 

3.3. Without relativity 

We now consider the evolution of Mercury's eccentricity without 
general relativity, i.e., described by the Hamiltonian H^l^ (24). 



Figure 4 shows the level curves of H^°] and H\f'', . The two 

& Mplan N,plm 

figures are similar to the subpanel a 20 of Fig. 1 within the region 
of low eccentricities {e\ < 0.6). But for higher eccentricities, two 
new fixed points appear in the Newtonian case with n = 20. It is 
then necessary to extend the expansion up to the order n w 50 to 
get the asymptotic topology of the phase space. Once the limit 
is closely reached, the result matches the relativistic ones well. 
The only difference is in the value of 6 g : for a given frequency 



/(SO) 




time (Myr) 

Fig. 5. Top: evolution of the fixed points of 7/ (50) without rela- 
tivity. Bottom: evolution of the eccentricity with time, using the 
simple model without relativity. 



gi, it should be increased by g r = 0.41"/yr with respect to the 
relativistic case (see Table C.l). 

Figure 5 shows the positions of the fixed points and the tem- 
poral evolution of the eccentricity driven by #Jyp] an - There is no 
significant difference with respect to the relativistic case. 



4. Spatial model 

The evolutions studied in the previous sections are perfectly reg- 
ular. Indeed, the Hamiltonian contains only one degree of free- 
dom and is thus integrable. The goal was to show that it is pos- 
sible to explain the large increase in Mercury's eccentricity ob- 
served in numerical simulations (e.g., Laskar 2008). Our aim is 
now to introduce an additional degree of freedom to make the 
system non integrable and to reproduce a chaotic evolution ex- 
pected in the vicinity of the separatrix of the g\ - gs resonance. 

To do so, we add inclination in our system, and focus on the 
term related to (gi - g$) - (s\ - S2). This resonant term is at 
present in libration, and was identified as one of the main source 
of chaos in the Solar System (Laskar 1990). It was indeed also 
computed analytically in (Laskar 1984), where it was identified 
as the major obstacle for the convergence of the secular pertur- 
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bation series. The Hamiltonian of the inclined problem is (see 
Appendix D) 



Hra 



1 



l-e\ 



pi\mo \a p 



4 

xr!; m {ei,A p , h , ip, m - ^5 > Q i - Q * ) ) . 



(29) 



where i p and Q. p are the inclination and the longitude of the as- 
cending node of the planet p, and I p and O* are the amplitude 
and the argument of the term precessing at the frequency S2 in 
the quasiperiodic decomposition of sin(/ ; ,/2) exp^Q^), respec- 
tively. In (29), both m* = #5/ + if* and Q* = s 2 t + <p* are func- 
tions of time t . To reduce the number of degrees of freedom from 
2.5 to 2, we apply transformations similar to those of the planar 
case (Sect. 2). We first add to the Hamiltonian the momentum 
T conjugated to the time t. The old variables of the subsequent 
canonical transformation are 



7=1- Jl-e 



-TJT\ , 



J = -yjl - ej(l - cos /1) , \p = -Qi 

f, t, 



(30) 



and the new variables are denoted (7, Am), (J, AQ), (T , t). We set 
Am = -6 - g 5 t - ip* = mi - m* , 

AQ. = -ijr- s 2 t - 0* = Qi - Q* , (31) 
t = t . 

The associated conjugated momenta are given by 
/ = 

/=-/, (32) 
t = -g 5 I -s 2 J + f . 



E 



0.5 



-0.5 




ex cos Aro 




150 200 250 
time (Myr) 



10(1 



Fig. 6. Top: trajectory of Mercury's eccentricity in the plane 
(e\ cos Am, e\ sinAtzr) for the inclined system. Bottom: evolu- 
tion of Mercury's eccentricity as a function of time in the same 
system. 



Since the new Hamiltonian H R ; nc is independent of t, its conju- 
gate momentum T is an integral of the motion. From the expres- 
sions of 7 and J (30), we thus get, up to a constant, 



g 5 + 2s 2 sin' 



l-el 



1 



7 



(33) 



To obtain the final Hamiltonian 77« ; nc in the spatial case, we add 
the control term and get 



1 2 

77R,inc - 77 Ri j nc + 2^g e i ■ 



(34) 



The Hamiltonian (34) now has two degrees of freedom. The vari- 
ables are (e\,Am) and (;'i,AQ). To perform numerical integra- 
tions, we use nonsingular rectangular coordinates 



k — e\ cos Am , 



h — <?i sin Am 



q — sin — cos AQ , p = sin — sin AQ 
H 2 1 2 



(35) 



With x - ~ e \> tne equations of motion read as (e.g. 
Bretagnon 1974) 



dk 
dt 

dh 
dt 



<977j nc 
+X ^ + 2x 



h I <977 inc 
1^— 



+ P 



-x- 



Bk 



k ( <977 ir 



+ P - 



dp 
d77 in , 



dp 



(36) 



dq 


1 dHinc 




L 9H mc 




dt 


4x dp 


~i 


( h dk 


dh 


dp 


1 8H mc 




(JH i0C 




dt 






( h dk 


dh 



Since Venus is the only planet with a significative ampli- 
tude associated to the frequency s 2 , it is also the only planet for 
which the inclination is taken into account. The initial inclina- 
tion of Mercury is taken from (Laskar 1990) where all the terms 
in its frequency decomposition, except the constant one, have 
been added. Doing so, Mercury's initial inclination is measured 
with respect to the invariant plane. 
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Figure 6 shows the results of a numerical integration of 
the inclined system. A large chaotic zone appears in the vicin- 
ity of the separatrix, as expected. Mercury's eccentricity now 
switches randomly between low-amplitude circulations and 
high-amplitude excursions, reaching values that are close to 1 . 

It should ne noted that the Solar System is not presently in 
this chaotic zone. As described elsewhere (Laskar 1990, 2008; 
Batygin & Laughlin 2008; Laskar & Gastineau 2009; Lithwick 
& Wu 201 1), the system is at present in a secular resonance (with 
resonant argument (gi — gs) — (si - S2)), in a state of slow chaos, 
with slow chaotic diffusion. This diffusion will quasi-randomly 
change the value of g\, which can then approach resonance 
with gs- The phase diagram of Fig. 6 describes the behavior of 
Mercury's eccentricity once this slow diffusion has brought the 
system into the vicinity of the g\ - gs resonance. 



5. Conclusion 

In this study, we developed a simple model to account for the 
increase in Mercury's eccentricity due to the resonance g\ = #5. 
This coplanar model is based on the expansion of the perturbing 
function with respect to the semimajor axis ratios, and exact in 
eccentricity. In the resulting Hamiltonian, we kept only one term 
in the quasiperiodic evolution of the outer planets' eccentricity. 
We found that this approximation is sufficient to reproduce an 
increase in Mercury's eccentricity up to 0.8. But we also noticed 
that it is necessary to extend the expansion up to the order n = 20 
and « * 50 in the relativistic and in the Newtonian cases, respec- 
tively. The explicit form of this secular resonant Hamiltonian is 
provided in Appendix B. 

The asymptotic topologies of the phase space are very simi- 
lar no matter whether the relativity is taken into account or not, 
so these two cases just depend on the value of the resonant fre- 
quency (Table C.l). In both cases, the system is a one degree of 
freedom system that is integrable with a separatrix (Figs, l.d, 4). 
The eccentricity of the trajectories in the vicinity of this sep- 
aratrix can rise to very high values, up to 0.8. Moreover, the 
timescale of Mercury's evolution is in very good agreement with 
the one observed in the numerical integration of the full model 
(Laskar 2008) (Fig. 3). With this integrable model, the behavior 
of the numerical solutions computed in (Laskar 2008; Batygin 
& Laughlin 2008; Laskar & Gastineau 2009) are understood, 
but not the transition from a low-eccentricity regime to a high- 
eccentricity regime. To obtain such a transition, it is necessary 
to include an additional degree of freedom in the system, which 
transforms the separatrix into a chaotic zone. 

Since we know that the resonant term associated with (gi - 
gi) ~ ( s i ~ s i) has large amplitude (Laskar 1984, 1990; Lithwick 
& Wu 2011), we added this single term in the spatial problem 
which corresponds to our second model (Sec. 4). As expected, 
in this nonintegrable problem, the separatrix is replaced by a sig- 
nificant chaotic zone where transitions from low-eccentricity to 
high-eccentricity regimes occur quasi-randomly, as observed in 
the full system, with maximal eccentricity close to 1 (Fig. 6). 
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Appendix A: Explicit expression of the Hamiltonian 
development in the planar case 

The secular Hamiltonian H v \ w of a restricted two-planet system 
where the massless body is on the inner orbit reads as 



A, lan = Yj ^r^°\eue p , m - m p ) 



(A.l) 



*P n=2 



The indices 1 and p refer to the inner massless body and to the 
outer massive planet, respectively; m p is the mass of the outer 
planet; a = a\ja p is the semimajor axis ratio; (ek)k=i,p and 
(•njk)k=i,p m ' e me eccentricities and the longitudes of periastron 
of the two planets. 

In (A.l), the T^ m {e u e p ,Tu) are given by (see Laskar & 
Boue 2010) 

<FZ'>(ei,e p ,'aT) = e n f n n_X (e x )X Q y (e p ) 

K«-D/2] (A.2) 
+ J] 2f lhq x;'"- 2 \eOX in+l) -"- 2 \e p ) cos((« - 2q)m) , 

where e„ — 1 if n is even and 0, otherwise 
(2q)\(2n - 2q)\ 



fn,q — 



2 2 "(q\) 2 ((n - q)\) 2 ' 



(A.3) 



and Xg' m (e) are Hansen coefficients. 

Now, we assume that the eccentricity of the outer planet re- 
mains low, and we only keep the linear terms in e p . Since the 

lowest power in eccentricity of X (n+l) '"'(e) is e'" 1 ', all the terms 
in the sum (A.2) are dropped, except those for which n-2q < 1. 
Furthermore, from the explicit expressions of the Hansen coeffi- 
cients (Laskar & Boue 2010), 



v—njn ^ 

~ (1 _g2)n-3/2 

[(n-2-m)/2] 



I 



(n-2) 



(\{m + {)\(n - 2 - (. 



! / e \ m+21 

-(m+2€))\ \2/ 



(A.4) 



for n > 2, one gets 



1 + 0(e 2 ) , 



(A.5) 



Then, the expression of T% simplifies as 



r^ 0) (eue p ,m) = e n f„ A X^(e l ) 

+(« - 1)(1 - e n )f n ^xf(er)e p cos m (A.6) 
+0(e 2 p ) . 

Substituting this expression into (A.l), one obtains 
Gm„ 



^plan — 



P ff (a, e\ ) - eie ; ,Q ff (o', e 2 ) cos^ - m p ) 



{A.l) 



+0(e 2 p ) , 



where 



P ff (a, e 2 ) = £ ff 2 V2 9 X*'V), 



(A.8) 
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and 



1 °° 

Q a (a, e 2 ) = --J] 2qa 2 ^f lq+u xlr l \e) 



Table B.l. Coefficients of the polynomials P(x) = £ p(X e (B.4) 
and Q(x) = £ q^x (B.5) computed up to a 50 . 



(A.9) 



9=1 



To derive explicit expressions of these two quantities, we use the 
analytical formulae of the coefficients X'^ m (e) for n > 2 (Laskar 
& Boue 2010), 



Xj"(e) = (-1> 



An + 1 — m) ! 



(«+ 1)! 

[(«+l-m)/2] 

Z 



(n + 1 — m)! 



n{m + €)\{n + l-m-2€)\ 



1+21' 



(A. 10) 



Moreover, we invert the sums and separate the cases I — and 
£ > 1 in P t /a, x). This gives 



Pi£r,x) = ]T/ 2M a 2? 

9=1 

°° ' °° (2^ +1)! 





p f x 10 s 


q t x 10 s 


I 


p t x 10 8 


q t X 10 8 


U 


j jOjo jyy 


IZlj 1UU 


1 J 


/OJL / jy 


44UZ oyy 


1 


C.Z 1 1 Q 7QH 

j 1 1 j / o / 


Z4j0 300 


1 A 
14 




Zl 1 / / oo 


Z 


1 n 1 OA AT7 
IV lOO O/ / 


ZjOo UJo 


1 J 


IVoj / /4 


/ / j / /z 


■2 
j 


1ZoV4 yjZ 


ZV03 /UU 


1 A 
10 


0/4 yt t 


O 1 1 AHA 

Z 1 1 oUo 


,-1 
4 


1Uo4j /ol 


IS. A A no /I 

J j44 UV4 


1 / 


1 71 TOO 

i IS ZSi 


41 Vj4 


j 


1 0?OR 7^fS 

1 UZ.UO / JU 


4^77 4^S 


18 


32433 


J o / J 


6 


10268 763 


5332153 


19 


4 303 


560 


7 


10780251 


6598 519 


20 


390 


35 


8 


11612770 


8057 388 


21 


23 


1 


9 


12557438 


9404 855 


22 


0.8 


0.03 


10 


13156622 


10023 739 


23 


0.02 


0.0002 


11 


12724433 


9279738 


24 


0.000 1 


0.0000006 


12 


10785 258 


7137 951 


25 


0.0000003 





(A. 11) 



'm\(2q+l -2€)\ 



v 2 'l 



(i)' ■ 



and 



We now define the polynomials P and Q such that 

#A<,plan = ^W.plan/Ai 

= -g 5 P (el) + g 5 Q (efje! cos(>i - m*) , 



(B.3) 



Q<la, x) = ^ ^ / 29+ i j? 

£=0 V ?=f 



(A. 12) 



where Ai = y]Gmoa\. The new polynomials P and Q are derived 
from P„ and Q„ through 



q(2q + 3)(2g + 1)! , 
£!(^ + 1)!(2^ + 1 -2^)! a 



It should be noted that in (A. 11) and (A. 12), the upper limits of I 
and g are infinite because we consider here the infinite expansion 
of the perturbing function in series of a. Nevertheless, when the 
Hamiltonian is truncated at the order a", the sums become finite. 
In P ff , the maximum value taken by I and q is [n/2], and in Q a , 
it is [(n- l)/2]. 

Appendix B: Numerical coefficients of the 
polynomials P and Q 

Here, we consider the system composed of Mercury perturbed 
by the seven outer planets from Venus to Neptune (p = 2, 8). 
The secular Hamiltonian of this problem, truncated at the first 
order in the outer eccentricities reads as 



g5 j-j mo a p \a p ) 



and 



Q(*) = — X — — Qo — ,x \xA p . 
85 mo a p \a p ) 



(B.4) 



(B.5) 



Their numerical values, summarized in Table B.l, have been 
computed up to the order a 50 from (Laskar 1990). 

Appendix C: Frequencies at zero eccentricity and 
corrections 

The precession frequency g\ of Mercury's perihelia computed 
either from the Newtonian Hamiltonian (20), 



Af,plan 



■z 



Gm, 



P =2 Q P 



a p 



(B.l) 



— .ejjcosfcji - vr p ) 

Then, using the quasiperiodic expansion of the variables 
(e / ,e' Er '') p= 2,8 and keeping only the terms at the fundamental fre- 
quency g 5 (A p e lvx i) p=2> 8, one gets 



#A<,pian = ~g5 yjl-e* ~ 8sP (el) + g 5 eiQ (efj cos Am , (C. 1) 
or from the relativistic Hamiltonian (23) 

g 5 ( V 1 ~ e \ + p ( e i) - e iQ ( e ?) cos Am ) 



H, 



(C2) 



Sr 1 



W.plan 



z 

p=2 



Gm, 



P a \^,el 



(B.2) 



-e\A p Qj^-,el Icosi cTj, - CT * ) 



is far enough from g$ a 4.249"/yr for the system not to be in 
secular resonance (see Table C.l). Nevertheless, due to the slow 
diffusion of the inner planets of the Solar System, this frequency 
is subject to small variations, and Mercury can eventually reach 
the resonance. To model this change in frequency, we include an 
additional term in the Hamiltonians: 6 g e\l2. 
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Table CI. Frequencies at zero eccentricity and corrections for 
different orders of expansion of the perturbing function. 



According to Laskar & Boue (2010, Eq. (B.30)), the term in 
factor of a" in the perturbing function reads as 



Newtonian case 



Relativistic case 



n 


Si f'Vvr) 


Sp C'/yr) 


pi CVvr) 


C'/yr) 


z 


i n 1 i 1 o 
j.y 1 1 lo 


U.JJ /04 


4.JZZ0J 


n nT/im 
-U.U/4U1 


5 


j.y i 1 lo 


U.JJ /04 


a moon 
4.JZZ0J 


n n~7/im 
-U.U/4U1 


A 

4 


4.V43oU 


n AO/i no 

-u.oy4yy 


^ T^^/IA 

J.JJJ4D 


1 1 nAA/i 
-1. 1U004 


c 

J 


4.y4joU 


n AQ/1QQ 

-u.oy4yy 


J. JJJ40 


1 1 nAA/i 
- 1 . 1UD04 


A 



C TOO AO 

J.jZoUo 


-l.U /yzO 


j. / jy / j 


1 zionoo 
- 1 .4yuyz 


1 


j.JZoUo 


i mniA 
-l.U /yzo 


j. / jy / j 


1 /innm 

-i.4yuyz 


Q 
O 


j.404oj 


i o i Am 
-I.zIoUj 


j.o /OjU 


1 AT7AO 

-i.oz /oy 


Cl 

y 


j.404oj 


i o i Ann 
-I.zIoUj 


J.O /OjU 


1 AnAn 

-i.oz /oy 


1U 


^ ^ 1 inn 
J.J IzUU 


1 OAT 1 O 
-l.ZOjlO 


j.yzjoj 


1 AT /lOQ 
-1.0 /40J 


i i 
1 1 


c c i inn 
j.j IZUU 


1 T AT 1 

-l.ZuJilo 


J.yZjOJ 


-1.0 /40J 


1Z 


Z CT70Q 
J.JZ/OO 


i nnm 

-i.z /yu / 


^ nine/I 
J.yjVj4 


1 Anmo 
-i.oyu /z 


i n 
1 J 


J.JZ /oo 


i nnm 

-i.z /yu / 


J.yjVj4 


1 Anmo 

-i.oyu 11 


1 A 

14 


C 1 /I 

j.j jj!4 


1 1 O/l T2 

-1.zo4jj 


^ ri/i /i on 
j.y44oU 


1 Ao^nc 

-i.oyjyo 


ID 


j.j jj!4 


1 1 O/l T2 

-1.zo4jj 


Z CiA a on 

j.y44oU 


1 Ao^no 

-i.oyjyo 


1 A 
10 


C CQ/1 QA 

J.J j4o0 


i 10 An/1 
-1 .Z50U4 


c Q/1 A 1 ^ 1 

J.y40jl 


1 AQ77n 

-i.oy / /u 


1 7 

1 / 


C CQ/1 OA 

J.J j4o0 


i 10 An/1 
-1 .Z50U4 


c Q/1 A 1 ^ 1 
J.V40J1 


1 AQ77H 

-i.oy / /u 


18 


5.53541 


-1.28659 


5.94706 


-1.69825 


19 


5.53541 


-1.28659 


5.94706 


-1.69825 




J.JJJJ7 


- 1 .Z.OU / / 




1 .U70H-Z. 


21 


5.53559 


-1.28677 


5.94724 


-1.69842 


22 


5.53564 


-1.28683 


5.94730 


-1.69848 


23 


5.53564 


-1.28683 


5.94730 


-1.69848 


24 


5.53566 


-1.28684 


5.94731 


-1.69850 


25 


5.53566 


-1.28684 


5.94731 


-1.69850 


26 


5.53567 


-1.28685 


5.94732 


-1.69850 


50 


5.53567 


-1.28685 


5.94732 


-1.69851 



Notes: g[ is Mercury's precession frequency computed at ej =0 with- 
out the correction S g e 2 J2 in each Hamiltonians. 6y represents the cor- 
rection S g that puts the system in resonance (gj = g 5 ). 



The values of 6 g putting the system in exact resonance at 
zero eccentricity are shown in Table C.l. One can observe that 
the correction is higher (in absolute value) when the relativistic 
precession is taken into account, which explains why relativity 
stabilizes the Solar System. 

Appendix D: Explicit expression of the Hamiltonian 
development in the spatial case 

Here we develop the secular Hamiltonian of an inclined re- 
stricted two-planet system following Laskar & Boue (2010). 
We use the same kind of approximation as for the planar case 
by considering only the linear dependency on the eccentricities 
and inclinations of the outer planets. We note ;'i and i p the in- 
clinations of the massless planet and of the perturber respec- 
tively, and Qi and Q. p are their longitudes of ascending node. 
We also note c\ = cos (;'i /2), si = sin (i\ /2), c p = cos (i p /2), and 
s p = sin (i p /2), from which we define 



/"* = [ c lCpS 2 + SxSpS 2 I 

v* = lcii p e ^ - iic p e ^ I 



rr } = y v v^-^x-^-^ep) 



j=0 q=0 



(D.2) 



xe i(«-2.s)n>] e \(n-2q)w p j 



Since we only consider the harmonics (wi - m p ) and ({tu\ - 
vj p ) - (Qj - ; ,)), we keep the terms such that q = n — s, and 
drop the others. Consequently, one sum disappears from (D.2), 
and it remains 



(D.3) 



rr - 1 (G ( ! l.^.,v t )x-"- 2 '( ei )x ( " +1) - 2 ' ! ( e;) ) 

xe i(«-2.s)Aw j 

Then, we extract the linear terms in e p . Since XZ (n+l) ' m (e p ) oc ej^ 1 ', 
we consider only the values of s such that n — 2s = 0, ±1. Using 
the asymptotic expressions of the Hansen coefficients (A. 5), we 
get 



^T = e^v*,v^>i), 

and 



(D.4) 



(0,0) _ fi{2n+l) flt ,^y2»tl,l L U< , J&w ^ (D 5) 



where cc means complex conjugate. Then, we use the definition 
of the Q functions (Laskar & Boue 2010, Eq. (B.31)). For s + q < 
n, one has 



with (Laskar & Boue 2010, eq. (B.32)), 

.00 1 (2s)\(2n-2q)\ 

q - s ,n- q -A > 2 1 " s\(n- s)\q\(n-q)\ 

x2c-i> 



(D.6) 



(2n-2s + k)\ 



(D.7) 



/.=o 



(2s-k)\(2n-2q-2s + k)\ k\ 



We thus have 



1 (2» )!(2»)! ^ ■ (2n + k)\ , (D ' 8) 
n! )4 Z.^ ^ (2 „ ■ ■ ■ - I' I 



2 4 " (n!> 
or in a more condensed form, 
Q<%>Qi*, v») = / 2n ,„F(2n + 1, -2n, 1; |v.D 



(2n-k)\(k\) 2 



(D.9) 



where f is the hypergeometric function and f pxi is given in (A. 3). 
In the same way, 



QTrlY = fi*f2n+uF(2n + 3, -2», 1; |v.|) 



(D.10) 



Now, we substitute the expressions of fi* and v, (D.l) into those 
of Ql n n (D.9) and of Q^+\ 10 )- Kee P in S onl y the linear terms 
in inclinations, 



ji t « c 2 e lAn + 2c 1 5 1 5 p 



(D.l) and 



|v»|* ~ sf - 2kcisf- l s p cos AQ , 



(D.ll) 
(D.12) 
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with AO = Qi - Q p , and retaining only the terms involved in the 
secular resonances, we get 

T£ 0) = f 2 n, n F(2n + 1, -2n, 1; $X^\e{) ; 

= 2f 2 „ +u {c\F(2n + 3, -2n, 1; sf) cos Am (D .13) 

+C! si s p G2„+\,„(s 2 l ) cos(Anr - AO.) j xf" +u \ei) ne p , 



where 

G2n+i,n(x) — 2F(2n + 3, — 2n, 1; x) 

-(1 -x)F'(2n + 3,-2n, l;x), 

and 



(D.14) 



ab 

F'(a,b,c;x) = —F(a+ \,b + l,c + 1; x) . (D.15) 

c 

The Eqs. (D.13) generalize the expression (4) obtained in the 
coplanar problem. The Hamiltonian of the inclined system is 
then 



Gm p 



a n r^ m (e u e p , i u i p , m - vr p ) , (D. 16) 



a P 

1 n=2 



with^, (0 ' 0) given in (D.13). 
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